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Abstract 

The vulnerability of an isolated network to cascades is fundamentally affected by its 
interactions with other networks. Motivated by failures cascading among electrical grids, 
we study the Bak-Tang-Wiesenfeld sandpile model on two sparsely-coupled random regular 
graphs. By approximating avalanches (cascades) as a multi-type branching process and 
using a generalization of Lagrange's expansion to multiple variables, we calculate the dis- 
tribution of avalanche sizes within each network. Due to coupling, large avalanches in the 
individual networks are mitigated — in contrast to the conclusion for a simpler model [36] . 
Yet when compared to uncoupled networks, interdependent networks more frequently suffer 
avalanches that are large in both networks. Thus sparse connections between networks sta- 
bilize them individually but destabilize them jointly, as coupling introduces reservoirs for 
extra load yet also inflicts new stresses. These results suggest that in practice, to greedily 
mitigate large avalanches in one network, add connections between networks; conversely, to 
mitigate avalanches that are large in both networks, remove connections between networks. 
We also show that when only one network receives load, the largest avalanches in the second 
network increase in size and in frequency, an effect that is amplified with increased coupling 
between networks and with increased disparity in total capacity. Our framework is appli- 
cable to modular networks as well as to interacting networks and provides building blocks 
for better prediction of cascading processes on networks in general. 

Keywords: sandpile models, random graphs, modular networks, branching processes, 
avalanches, cascades, self-organized criticality. 
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1 Introduction 



The interdependence among systems is paramount. A system that is stable in isolation, for 
instance, may lose stability when coupled to another system. By contrast, some systems are 
largely useless unless coupled to others; important examples of this occur in biology, in which 
modular components — organs, tissues, cells and organelles — function properly only in churning, 
interdependent synchrony. Such interdependence among systems can be detrimental, however, 
as small failures can escalate to catastrophe. Just one failed element within one module of a 
system — a malignant lymph node in a human body, an eradicated autotroph in an ecosystem, 
a downed power line in an electrical grid — can spread to other connected subsystems and, via 
feedback loops and percolation on a complex web of links, cascade to system-wide failure. The 
connections among systems can enhance or inhibit, hasten or delay such cascades. 

One of the best examples of coupled systems prone to cascading failures is infrastructure [Ml 
Il4" ll42[|4"l"] . There is no denying the increasing interdependence of modern infrastructure: water, 
gas, Internet packets, financial transactions, phone calls and electrons move on regional, national 
and global networks, responding to demand and depending on one another for proper function. 
When an element of one of these networks is overwhelmed, its neighboring elements may pick 
up the slack — they are often engineered to do so — but when they cannot, the failure can spread, 
both within that network and to others. 

Striking examples of failures cascading among interdependent infrastructure abound. In 
1998, for instance, a telecommunications satellite over the western United States malfunctioned, 
crippling the communication network, which in turn disrupted the transportation network be- 
cause gas stations could not process credit card transactions and airports lacked precise weather 
information [UJ. In Italy in 2004, the failure of a telecom node paralyzed nearly all national 
telecommunication services, in turn halting most Italian financial transactions, postal deliver- 
ies and flights [UJ. More recently, the ash-covered skies over Europe in the aftermath of the 
eruption of the Eyjafjallajokull volcano in Iceland on April 14, 2010, halted nearly all flights 
in western Europe, which not only disrupted cultural, sporting, military and diplomatic events 
but also affected economies across the globe due to their unprecedented interdependence. It is 
remarkable the extent to which failures and other information ripple among today's increasingly 
intertwined economies, penetrating the far corners of the globe: thousands of flower farmers in 
Kenya, for instance, were laid off because their harvests could not be shipped by air to the UK 
in the aftermath of Eyjafjallajokull, layoffs that in turn affected the Kenyan economy [23j. 

Such interdependence motivates analytic models of coupled networks. In these models, 
each infrastructure system is a network, yet nodes in one network may depend on nodes in 
other networks. For example, the World Wide Web is a virtual network that runs on the 
physical Internet of routers, which in turn depends on the electrical grid for power, all three 
of which are (increasingly) crucial to the finance, transportation, health and other systems. 
The failure of nodes in one network may cause nodes in other networks to fail, which in turn 
impairs nodes in other networks, and so on. This recursive structure of the failures hints at 
the possibility of rigorous, mathematical analysis (e.g., [12J ) . A better understanding of the 
characteristics of coupled networks that facilitate or inhibit spreading processes on them may 
help to answer the difficult question of design: How do we fortify interdependent infrastructure 
against catastrophe? Furthermore, how can we effectively spread desirable "pathogens", such 
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as information and ideas, among interdependent networks? 

The idea of interacting networks is closely related to the concept of modules or communities 
within networks. In the classic definition, modules are subgraphs with significantly more edges 
than expected within the subgraph and few edges to the rest of the network |19j . Alternate 
definitions of modules have been studied that capture, for example, hierarchical structure [13] 
or similarity in patterns of connectivity [38]. For a recent comprehensive review of community 
structure, see Ref. [TS]. Despite the extent of research on finding modules in networks, their use 
in applications so far has been limited. Our aim is to use knowledge of the modular structure of 
isolated networks (or of the structure of interactions among distinct networks) to better predict 
the behavior of the whole system. 

The "electrical grid" of the United States, for example, is not a grid but instead consists of a 
collection of over 3,200 distinct local and regional electric utilities [47J. A small number of these 
are federal utilities, while the vast majority are either investor-owned or publicly-owned. Each 
individual utility is an independent entity, with myriad economic and physical considerations 
dictating its structure and operating policies, yet the entities connect together to form the 
"grid". Historically the grid has comprised a collection of loosely interconnected, local systems, 
but the level of connectedness, length scales of interactions, and the number of small, distributed 
power sources are increasing [5], Given this increasing level of diversity and interconnectivity, 
it is important to understand at a fundamental level the impact of interactions on cascading 
failures. 

With coupled electrical grids in mind — though not limited to this application — we study the 
Bak-Per-Wiesenfeld (BTW) sandpile model [7] on interacting networks. In the sandpile model, 
grains of sand are randomly dropped on nodes, each of which possess an innate capacity to hold 
sand. Whenever a node's load exceeds its capacity, it "sheds" its load to its neighbors, and 
any nodes now exceeding their capacity shed their load synchronously at the next time step. 
This process of unstable nodes toppling continues until equilibrium is restored. In this way, 
dropping a single grain of sand can cause an avalanche, often small but sometimes spanning 
the entire network. Sandpile models on networks are toy models that do not perfectly capture 
the structure and dynamics of electrical grids, but they share enough properties to be relevant 
by elucidating the space of possible behaviors in simple models. Moreover, sandpile models are 
simple enough to find applications in other systems that bear and shed load and that yield 
cascades with the long correlations and power-laws that characterize self-organized criticality 
(e.g., neuronal networks [HI [28]; forest fires, earthquakes, landslides, financial markets [6]). 

There is a large body of research on cascading failures in infrastructure (e.g., [34 4 1441 l36"l \A2 \ 
[T5l |4"TI 145]). but only recently has the effect of the interdependence among systems garnered 
considerable attention. Newman et al. have studied two such models: CASCADE, a sandpile 
model in which load (sand) is shed uniformly to nodes in two networks in a non-local way |36j ; 
and DSCM, a probabilistic model of random failures and repairs in which neighbors of failed 
nodes are more likely to fail [36] ■ The CASCADE model, approximated by a branching process, 
demonstrates how the coupling of networks shifts the critical points of global cascades. However, 
to reach this conclusion analytically they make the unrealistic (but mathematically convenient) 
assumption that load is shed globally to every node rather than locally, so that critical points 
are simply maximal eigenvalues of the network coupling matrix. Meanwhile, DSCM (short 
for "coupled dynamic complex system model") is a simple model that exhibits self-organized 
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criticality — i.e., the system dynamically arranges itself at the critical point, where it exhibits 
long time correlations and power-law distributions of failure size. Newman et al. conclude that 
introducing coupling between systems makes them more vulnerable, though the models are 
simplified and deserve more sophisticated, realistic studies. 

A model due to Bulydrev et al. p2] is also simplified for the sake of analytical solutions. 
They consider two networks, each with the same number of nodes, and these networks inter- 
act in that each node in one network is connected uniformly at random to one node in the 
other network. Initially, a random subset of nodes in one network fail, and then a cascade of 
failures ensues, in which internal edges that connect disconnected components (i.e., maximally 
connected subgraphs) in the other network are deleted. Although the assumptions regarding 
the network coupling are somewhat unrealistic (every node has one neighbor in the other net- 
work, chosen uniformly at random) and although the failure mechanism is very specific (delete 
internal edges that connect disconnected components), the model nevertheless demonstrates 
the surprising catastrophes possible in interdependent networks. Moreover, considering the in- 
terdependence between networks can upend longstanding notions about isolated networks. For 
example, Buldyrev et al. found that for coupled networks with their failure mechanism, broader 
degree distributions are more vulnerable to random failure |12j . in contrast to the conclusion 
for isolated networks that broader degree distributions are more robust to random failure [2]. 

In this paper we analyze cascades in sandpiles on interacting networks by approximating 
them by multi-type branching processes. Here we consider locally tree-like networks, as are the 
case for random graphs generated using the configuration model as in Ref . [39] . Studies of the 
electrical grid suggest that they are tree-like (low clustering coefficient) and have narrow degree 
distribution: two small-scale studies of the topology of nodes constituting a small regional 
electric grid suggest that the connectivity of nodes obeys an exponential distribution [H Q], 
and a recent extensive study corroborates the exponential nature of the connectivity as well 
as a small non-zero clustering coefficient (C = 0.071) [26J. The framework we develop is 
for arbitrary degree distribution but requires tree-like graphs (with C ~ 0). Extensions to 
social and other networks that are not tree-like would build off recent work on random graphs 
with arbitrary subgraphs [Ml E3 EZ]- We solve the multi-type branching processes for the 
avalanche size distributions using an old result in the mathematics literature, a generalization 
of Lagrange's expansion due to I. J. Good |22| . as well as numerically simulate a situation not 
yet described by our analytic framework: asymmetric load applied to one network rather than 
to both. By studying processes on modular or interacting networks using multi-type branching 
processes and by using knowledge about the connection structure within and between modules 
(as in [201 HH1 HH1 OS HQl [3]), we begin to develop a more complete picture of the behaviors of 
complex, heterogeneous networks. 

We show that introducing coupling between networks stabilizes them individually, in the 
sense that large avalanches are mitigated and small avalanches are amplified, yet destabilizes 
them jointly, in that avalanches that are large in both networks become more likely — compared 
to the null hypothesis of two uncoupled graphs. The finding that connecting networks stabilizes 
them individually contrasts to the result in [36] for a simpler sandpile model, and it assuages 
the warnings in [T2] about the catastrophic cascades of failed connectivity in a model of coupled 
networks. On the other hand, we also find that introducing connections between networks can 
destabilize them individually in some circumstances: if only one network receives external load, 
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the other network suffers large avalanches that increase in severity and frequency with increased 
coupling between networks and with increased disparity in relative capacity. This suggests an 
arms race to increase capacity of a network to fortify against cascades inflicted by neighboring 
networks. 



1.1 Classic BTW sandpile models on lattices 

Introduced in the late 1980s [7J, the Bak-Tang-Wiesenfeld sandpile model is a well-studied toy 
model of cascades that exhibits self-organized criticality, power laws and universality classes. 
In a classic version of the model on a finite, two-dimensional lattice, grains of sand are dropped 
uniformly at random on nodes in the lattice, and whenever a node contains four or more grains 
of sand, it sheds one grain to each of its four neighbors at the next time step. The lattice has 
open boundaries, so that sand shed off the boundary of the lattice is lost, which prevents the 
system from becoming inundated with sand. Various measures of the size, area, and duration 
of avalanches follow power- laws, and these variables relate to one another via power laws |10j . 

A few variants of the sandpile model on lattices can be solved exactly if the shedding rules 
have abelian symmetry |10| . More recently, sandpile models have been studied on (isolated) 
networks, including Erdos-Renyi graphs [11, 33j, scale- free graphs [21], and graphs generated 
by the Watts-Strogatz model on one-dimensional [30] and on two-dimensional [H] lattices. A 
common question explored using asymptotic calculations and computer simulation is: Under 
what conditions on the network structure is the avalanche behavior "mean-field" — i.e., approx- 
imately that of a complete graph — which corresponds to the avalanche size distribution being 
a power law with exponent 3/2. 



1.2 BTW sandpile models on arbitrary networks 

The most natural choice for the "thresholds" or "capacities" for each node in a network is its 
degree, so that nodes shed one grain to each neighbor (or, more precisely, one grain along each 
outgoing edge, since it may have parallel edges or self- loops) [21]. Another choice for sand 
thresholds is uniform |21| 133]: for example, leaves (i.e., nodes with degree one) have threshold 
one, while all other nodes have threshold two. However, uniform threshold suffers from two 
drawbacks: shedding becomes ambiguous because nodes must shed to a random subset of their 
neighbors rather than one grain to every neighbor, and a large portion of the novelty of sandpiles 
on networks (rather than on lattices) is the disparity among nodes in the amount of sand shed. 
As a result, here we choose thresholds of nodes to be their degrees. 

To make this explicit, the sandpile dynamics on networks are as follows. With the graph 
fixed, we add grains of sand to nodes in the network chosen uniformly at random. Whenever 
the sand on a node exceeds its total degree, the node "topples" and sheds one grain to each of 
its neighbors. (More precisely, it sheds one grain along each edge, so if it has multiple edges to 



a neighbor, then it sends as many grains to that node, but as discussed below in Section 2.2 
parallel edges in the sparse random graphs considered here are rare.) If any new nodes become 
unstable, they all shed synchronously at the next time step. (If a node has strictly more grains 
than its degree, it still sheds one grain to each neighbor, and the leftover grains remain on the 
node; in our studies, the chance that a node has at least twice as many grains as its degree — and 
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hence has to shed multiple grains per neighbor — is negligibly small.) Topplings continue until 
no node exceeds its capacity — i.e., equilibrium is restored — whereupon the process repeats. 

To ensure that the network doesn't become overloaded with sand, whenever a node sheds 
its sand, each grain is deleted independently with probability /, which we call the dissipation 
rate of sand. The dissipation of sand in sandpile models on networks is the analogy of open 
boundary conditions in classic BTW sandpile models on finite lattices, in which grains shed off 
the boundary are lost. This parameter can profoundly affect the duration and size of avalanches. 
In [21], they used / = 0.0001 for a scale- free network with 10 4 nodes, but they did not mention 
the sensitivity of the avalanche size on /: decreasing / prolongs avalanches, which can add a 
"hump" to the avalanche size distributions at large avalanches. We discuss how to choose / in 
Section 13.41 below. 

2 Random graph model of two coupled networks 

Here we study the sandpile model on two interacting networks, labeled a and b, which have 
their own internal (or intra-)degree distribution and which are sparsely coupled by edges, called 
inter-edges. (Throughout we use the prefixes intra- and inter- to refer to edges and degrees 
within or between the two networks, respectively.) The connectivity of interacting networks was 
studied in [32]. Each node has an intra- and inter-degree, which are the numbers of neighbors 
within its network and in the other network; we choose the thresholds of sand to be the most 
natural one, the total degree. 

2.1 Inter- and intra-degree distributions 

Since degrees play a key role in avalanche dynamics, we study in detail the degree distributions 
of interacting networks. Each node has an integer number k a > many a-neighbors and kb > 
many 6-neighbors. Each network, a and 6, is characterized by a degree distribution on Z> : 

Pa(k a ,kb) = fraction of a-nodes with k a a-neighbors and kb 6-neighbors, (1) 
Pb(k a ,kb) = fraction of 6-nodes with k a a-neighbors and kb 6-neighbors. (2) 

The associated generating functions are 

oo 

G a (uj a ,uj b )= ^2 Pa(k a ,k b )uj^ a uj^ b , 

ka >k o = 

oo 

Gb(u a ,Ub)= ^2 Pb(k a ,k h )u% a U} b b , 

k a ,k b =0 

for td a ,Ub G C with \uj a \, \wb\ < 1 (i-e-> in the bidisc). 

Definition 1 We say that the marginal p aa {-) = X^=o Pa(' > ^b) * s the intra-degree distribution 
of network a, while the marginal p a b(') = =oPa{^a, •) is the inter-degree distribution of 
network a. The intra- and inter-degree distributions of b (pb a andpbb) are defined analogously. 
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It is convenient to assume that the intra- and inter-degrees of nodes are independent, but this 
is rarely the case in real- world networks; the connections between different infrastructures may 
occur most frequently between nodes of low degree, for example. 

Definition 2 The intra- and inter-degree distributions are independent if the degree distribu- 
tions can be written as a product of two probability distributions on Z>q: e.g., 

Pa(k a ,h) = Paa{k a )Pab(h), 
Pb(k a ,h) = Pba{ka)Pbb(h)- 

Remark 1 It may be that the degree distribution of one of the networks is independent while 
that of the other is not. 

2.2 Generating interacting networks from their degree distributions 

We generate coupled random graphs from their degree sequences using a simple extension of 
the configuration model. However, this necessarily changes the inter-degree distributions in a 
subtle way, due to the fact that the edges between networks are undirected. Below we make 
precise this effect of conditioning on the event that there must be equally many edges from a 
to b as from b to a, which holds in general for bipartite undirected random graphs and any 
interacting undirected graphs. 

The network generation works as follows. First, each node independently draws its pair of 
intra-degree and inter-degree, (k a ,kb), from its network's degree distribution, p a or pf,. Such 
i.i.d. degree sequences (one for each network) are drawn until the sum of the degrees within 
each network is even and the sum of the inter-degrees from a to b equals the sum of the degrees 
from 6 to a (since the edges are undirected). (Note that this second requirement conveniently 
vanishes for directed inter-edges.) Once valid degree sequences are drawn, each node's "half- 
edges" (or "edge stubs") of the two flavors (namely, "toward a" and "toward 6") are wired 
randomly, as in the configuration model (potentially with correlation between degrees, which 
we do not consider here). 

Three features deserve attention: parallel edges, few short cycles (i.e., it is locally tree-like), 
and the effective inter-degree distributions. First, the configuration model permits self-loops 
and parallel (multiple) undirected edges between two nodes. Parallel edges can be a desirable 
feature for a model of cascades because it allows variable amounts of connectivity: effectively 
they are single edges with integer weight. For sandpile avalanches, parallel edges mean a node 
could shed more than one grain to a neighbor. However, for the large, sparse random graphs 
considered here, parallel edges are rare. 

Second, we note that this configuration model generates locally tree-like graphs, which makes 
them amenable to branching process approximations. Many infrastructure networks are locally 
tree-like; electrical grids, for example, have few small loops, but they must have large loops 
since the electrons cannot disappear. By contrast, many social networks are emphatically not 
tree-like because of the small loops caused by transitivity: "friends of friends tend to be friends" . 
Such small loops and small cliques are exceedingly unlikely in tree-like graphs generated by the 
configuration model, since the graphs are sparse. Percolation and contact processes on clustered 
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networks |35[ [37] and ones with arbitrary distributions of small subgraphs [29] have recently 
been studied. 

Finally, we calculate the effective inter-degree distributions, which are not the input inter- 
degree distributions because the undirected edges between the two networks require that a has 
as many edge stubs toward b as b has toward a. Intuitively, if the inter-degree distributions 
p a b(-) and Pb a {-) have much overlap, then the effective inter-degree distributions hardly differ 
from the input ones. But if they profoundly differ — for example, two Poisson distributions 
with very different means — then the valid degree sequences, which condition on the number of 
inter-edges agreeing, have a very different distribution than the input degree distributions. For 
this example of two Poisson distributions with different means, the effective means would be 
much closer to one another, since valid degree sequences occur more frequently between the 
two means. 

Lemma 1 Let X,Y be the random variables for the inter-degree sequences of a and b, which 
are vectors of length N a ,Nb > 0, respectively. We denote Ex = Y7i=i x i t° be the sum of the 
entries in the vector. Then the effective inter- degree sequence for a is 

Pt(X = k\EY = Ek) = Pr{X = k) p*f b (Ek) (3) 
where p* b ^ b {■) is pb a convolved Nb many times. 

Proof. Using the independence of X and Y and by summing over all t G % Nb , we have 

Pv(X = k, EY = Ek) 



Pr(X = k | EY = Ek) 



Pr(EY = Ek) 
£fPr(X = k, Y = e, He = Efc) 

Pr(Ef = Ek) 
Pv(X = k) Pr(f = £, Yl= Ek \ X = k) 

Pr(Ey = Ek) 
Pr{X = k)Y^ Pr(T = t\X = k,EY = Ek) 



(4) 



Pl(X = %) (p ba *Pba * -*Pba) (Sfc) ■ (5) 



Nf, many convolutions 



Eq. Q follows from the independence of X and Y, and Eq. ^ follows from the definition of 
the convolution. □ 



Lemma [T] says that the effective inter-degree distribution of a is a product of the input 
inter-degree distribution p a b of a with the inter-degree distribution pb a of b convolved iV;, times 
and evaluated at the sum of the entries of the input to p a b- In the literature this effect is 
often overlooked in theoretical calculations. In practice, when generating bipartite graphs (or 
other interacting networks) from degree distributions using the configuration model, a common, 
quick solution is to draw degree sequences from their distributions, and then repeatedly choose 
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a node uniformly at random from anywhere in the network and re-draw its degree until the 
degree sequences are valid. However, this method suffers from the above effect, which often is 
subtle but can be substantial if the degree distributions have "little overlap" . The inter-degree 
distributions considered here have "much overlap", so the effective inter-degree distribution is 
approximately the input one, and the correction factor in Eq. §5§ can be neglected. 

3 Avalanche size 

Most grains of sand dropped onto the network do not topple any nodes; instead, they simply 
increase a node's load, but not beyond its capacity. Some grains of sand topple a node that in 
turn may topple a few others. Even fewer trigger large avalanches that topple nearly the entire 
network before equilibrium is restored. We are interested in the asymptotic distribution of the 
sizes of avalanches, after many grains of sand have been dropped. 

3.1 Two measures of avalanche size: topplings and sheddings 

To that end, we measure the size of avalanche in two ways: (1) the numbers t a ,tb of toppling 
events in a and in 6, respectively, and (2) the numbers of grains of sand that are shed from 
one network to itself or to the other network. For short we call the former "a-topplings" and 
"6-topplings" , and we call the latter "od-sheddings" , where o,d G {a, b} are the "origin" and 
"destination" networks of the shedded grain of sand. For example, the toppling of an a-node 
that has k a = 5 neighbors in a and kb = 1 neighbor in b counts as one a-toppling, five aa- 
sheddings and one a6-shedding. For a wide range of sandpile models, the various measures of 
avalanche size — size, area, perimeter, duration, maximal distance, radius of gyration, etc. — 
scale against each other in the form of power laws |10j . which suggests that it suffices to study 
just a couple of measures. 

Since the two networks have different degree distributions, the distributions of the number 
of topplings depends on the network in which the first grain is dropped. As an example, a grain 
dropped in a dense network a, weakly coupled to a sparse network b, would likely cause larger 
avalanches in b — compared to avalanches begun in b — since a has more total capacity, so its 
large avalanches — the most likely ones to cross the sparse inter-network edges — overwhelm the 
lower-capacity network b. Thus we have two distributions of toppling size, with the subscript 
indicating the network in which the avalanche begins: 



These toppling size distributions ('s' for "size") count the initial toppling, and so they satisfy 



since, for example, an a-node cannot topple 6-nodes unless it topples first. Note that s a and 
Sb are probability distributions in the asymptotic limit, in the sense that they are frequencies 



Sa(ta,tb) = prob. a grain dropped in a causes t a topplings in a and % in b, 
Sb{ta,tb) = prob. a grain dropped in b causes t a topplings in a and % in b. 



s a {0,t b ) = 
s b (t a ,0) = 



V t b > 1, 

V t a > 1, 
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of avalanche sizes after the network has undergone many cascades. For each avalanche size 
distribution, we define the associated generating functions 

oo 

S a {T a ,n)= 8 a (t a ,t b )T^T^ , 

t a ,t b =o 

oo 

S b {T a ,T b )= Sb(ia,t6)r* a T* b , 

t a ,t b =o 

for r a , neC with T a ,n < 1. 

For interacting networks, sheddings rather than topplings are a more natural way to math- 
ematically study avalanche size, for reasons discussed below. Recall the definition of an od- 
shedding: 

Definition 3 An o<i-shedding is the event that a grain of sand is shed from network o G {a, 6} 
to network d E {a, b}. (We use 'o' for "origin", 'd' for "destination".) 

In order to calculate the distributions s a ,Sb of the number of topplings in a and b, we must 
first calculate the distributions of the number of sheddings of the four origin-destination types, 
aa, ab, ba, bb. From these shedding distributions we then calculate the avalanche size distribu- 
tions S a ,Sb- 

We denote by p aa , p a b, PbaiPbb the probability distributions of the number of od-sheddings of 
the four types (aa, ab, ba, bb) in an avalanche caused by an initial shedding of the type specified 
in the subscript. Like the avalanche size distributions s a , Sb, these are defined asymptotically, 
after many grains of sand have been dropped on the network. The probability generating 
functions V d associated to the shedding distributions p a d (where o,d G {a, b}) are 

oo 

V d{o- a a,o- a b,a ba ,o-bb) = ^ Pod(raa,r a b,rba,rbb)o-l a a all b a r b b a a a2 b , 

r aa ,r ab ,r ba ,r bb =0 

for o- a a,o~ab,o~ba,o~bb £ C with absolute value < 1. 

3.2 Sheddings are topplings on the line graph 

Whereas cascades of topplings correspond to cascades of nodes in the original graph, cascades 
of sheddings, in a sense formalized below, correspond to cascades on the nodes of the line graph 
of the original graph. (Recall that the line graph [23] of a graph G is a graph with nodes 
labels given by edges in G, and two nodes in the line graph are connected if and only if the 
corresponding edges in G share a common vertex.) 

Now we make the connection between sheddings and the line graph precise. For a given 
avalanche let the toppling digraph be the directed graph in which there's an edge from node u 
to node v if u shed at least one grain of sand to v in the avalanche. Similarly, let the shedding 
digraph be the directed graph in which two nodes labeled '« — > v' and l w — > x' are connected 
if v = w and a shedded grain from u to v caused v = w to topple and to shed a grain of sand 
to x — i.e., if u toppled, which toppled v = w. 
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Proposition 1 If the graph is a tree, then the shedding digraph is the line graph of the toppling 
digraph. 

Proof. If node u topples, which eventually topples v, then every node on the unique path 
between u and v toppled in succession because the graph is a tree, and hence the successive 
directed edges along the path are connected in the shedding digraph. Conversely, if u is not 
connected to v in the toppling digraph, then — again since the graph is a tree — it never occurred 
that v toppled immediately after u toppled, and so for any node y there does not exist an edge 
from (u, v) to (v , y) in the shedding digraph. □ 



Remark 2 If the toppling digraph is only approximately tree-like, then this becomes less precise, 
but in a sense the shedding digraph is approximately the line graph of the toppling digraph. There 
may be a way to encode the avalanche in a weighted digraph that preserves all the information 
about "who toppled whom" in the avalanche. 

3.3 Branching process approximations 

Cascades in networks can be approximated by a multiplicative branching process if there are few 
loops of small or intermediate size. Said differently, a growing avalanche is approximately a tree 
with branches that grow and terminate independently. This approximation works well for locally 
tree-like networks, such as the random graphs generated here using the configuration model; 
since the edges are sparse and wired uniformly at random, the chance of closing a triangle (or any 
other small loop) is negligibly small. Two loosely coupled networks appear to be approximately 
tree-like as long as the coupling is sparse, as verified for coupled Erdos-Renyi and power-law 
graphs in [32J. For two coupled networks, the avalanches can be approximated by two-type and 
four-type Galton- Watson processes for the distributions of topplings and sheddings, respectively, 
which we formalize next. 

The fundamental tool in branching processes is the "children" or "branch" distribution, the 
distribution of the number of "children" objects of the various types generated by a "parent" 
object of a certain type. Many probabilistic properties of the tree generated by an initial 
population are calculated using the branch distribution (or even just its moments). It turns 
out that it is more natural to write down the branch distribution of sheddings rather than 
topplings. 

3.4 Shedding branch distributions q a a,Qab:QbajQbb 

To calculate the shedding distributions p aa , Pab, Pba, Pbb, we first must determine the distribu- 
tions q aa ,QabiQba,Qbb of the number of "children sheddings" caused by a "parent shedding" of 
the four types (aa, ab, ba, bb) specified in the subscript. 

Definition 4 The branch distributions q d( r da, r db), where od G {aa,ab,ba,bb}, are the proba- 
bilities that a grain of sand shedded from an "origin" node in network o G {a, b} to a "destina- 
tion" node in network d G {a, b} topples the destination node, which in turn sheds one grain to 
each of its r^ a many a-neighbors and r^b many b-neighbors. 



11 



Said differently, q d{r d a-,rdb) is the chance that an od-shedding causes r^a many da-sheddings 
and many dfr-sheddings. In the language of Galton- Watson processes, q d{ r da, r db) is the 
distribution of the number of "children sheddings" of types da, db G {aa, ab, ba, bb} caused by a 
"parent" shedding of type od G {aa,ab,ba,bb}. The "children" sheddings necessarily originate 
in network d G {a,b}. 

To determine the q Q d, consider a grain of sand that has just been shedded from an "origin 
node" to a "destination node". First we estimate the chance that this grain of sand topples 
the destination node (the node to which it is shed). Empirically, similar to |21j we find that 
the amounts of sand on nodes are approximately uniformly distributed from zero to one less 
their degree asymptotically — there is no typical amount of sand — so the probability that the 
destination node topples is the probability that it has one fewer grain of sand than its total 
degree: l/(k a + k\,)- (Later we mention correction terms that account for toppling due to 
receiving multiple grains of sand at once, which is possible only if there are loops in the network, 
and so they can be neglected for the locally tree- like networks considered here.) 

We are interested not only in whether the destination node topples but also in how many 
grains it sheds to each network, which is the destination node's intra- and inter-degree. The 
asymmetry of the coupling between networks implies that the degree distribution of the destina- 
tion node depends on the network to which the origin node belongs. For example, suppose the 
inter-edges only connect intra- hubs (the nodes with high intra-degree); then a grain traveling 
an inter-edge is more likely than one traveling an intra-edge to land on a hub (even though hubs 
are already likely to be found at the end of intra-edges!). This is why determining the branch 
distributions q a d of sheddings is easier than determining the branch distributions of topplings: 
we must carry two indices in the subscript to use the fact that the distribution of the degree of 
the destination node in network depends on the network to which the origin node belongs. 

If the networks a and b are connected by sparse edges without intra-degree correlations, 
then — to good approximation — a grain traveling from, say, network a to network b arrives at 
an inter-edge stub in b chosen uniformly at random. Similarly, if the intra-edge stubs are joined 
without degree correlations, then a grain traveling an intra-edge is equally likely to arrive at 
any other intra-edge stub. Since a d-node with r do > 1 many oneighbors is r do times more 
likely than a (f-node with one oneighbor to receive a grain of sand, for o,d G {a, b} we have 

, x rdoPd{rda,r d b) 1 f . ^ n 

q d{rda, rdb) = t, — r ; for r da + r d b > 0, (6) 

\kdo) r da + r db 

and we set 

<3W(0,0) := 1 - ^2 q d(rda,rdb)- (7) 

rda+rdb>0 

The denominator {k do ) := d Wo GdO-, 1) is the expected number of edges from d to o; it normalizes 
the numerator r<fo£><z(?"da> i~db) to be a uniform probability distribution over the edge stubs from 
network d to network o. Eq. Q is the probability that the destination node does not topple — 
i.e., that it has fewer grains than one less than its total degree. 

The justification of Eq. ^ is as follows. Define the following four random variables for a 
node in network b: 
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1. K a € Z>o, the node's a-degree 



2. Kf, £ Z>o, the node's 6-degree 

3. G £ Z>o, the number of grains on it 

4. i? £ {a, 6}, whether it receives a grain from network a or from b 

Next we define, for network b, the probability distribution Pi(K a , K b , G, R) on Z>o x Z>o x 
Z>o x {a, 6}. Then the chance that a 6-node receives a grain from an a-node and topples its 
grains to k a a-neighbors and kb 6-neighbors is 

Qab(k a , k b ) = Pr(K a = k a , K b = k b , G = k a + k b - 1, R = a) 

= Pr(K a = k a , K b = k b ) ■ Pr(G = k a + k b - 1 \ K a = k a ,K b = k b ) 
■ Pr(R = a\ K a = k a ,K b = k b , G = k a + k b - 1) 
1 K 



= Pb(k a ,k b ) 



k a + h YaLi T,T=o iPbihJ 



Remark 3 If the degree distributions of networks a and b are independent, meaning they can 
be written as products of intra- and inter- degrees 

Pa(k a ,k b ) = Paa(k a )p ab (kb), 
Pb(k a ,k b ) = Pba{ka)Pbb(h), 

then the branch distributions simplify to 

, s r do p da (r da )p db (r d b) 1 

qod{r da ,r d b) = T , — , , 

{kdo) r da + r db 

where (as above) we denote averages by (•) (e.g., (k aa ) = Yli^i Waa(i)))- As above, these 
equations hold for k a , kb not both zero; at the origin they equal one minus their values everywhere 
else in Z> . 



Note that the branch distributions q od are the distributions of sheddings at the next time 
step, so they depend on only two inputs — k da and k db — because a ci-toppling can only cause 
da- and (ifr-sheddings. By contrast, the shedding distributions p od count the total number of 
aa-, ah-, ba- and 66-sheddings in avalanches initiated by an o<i-shedding. Since an od-shedding 
could lead to sheddings of any type in the ensuing avalanche, the p od depend on four inputs: 

r aa, fab, r ba, r bb- 

If the inter-edges are assortative, disassortative, or some other pattern deviating from in- 
dependence, then the expressions for q ab and q ba need to change. For example, if the edges 
between the networks are assortative and the node that's shedding has high intra-degree and is 
in a, then q ab should be larger for nodes with high intra-degree. This needs a new, more general 
framework than the q od above, since we need to take into account the intra-degrees of the two 
nodes, which do not appear anywhere in our equations for the g's. Similarly, if a and b have 
intra-assortativity (or some other correlation), then we need to modify q aa and q bb to account 
for these correlations. 



13 



The q d are only approximations to the avalanche dynamics, and one could add correction 
terms for greater accuracy. First, if the network is not approximately tree-like (such as the 
Watts-Strogatz model on ID lattices or random graphs with triangles added), one could add a 
correction for loops by computing the chance that a node is two grains away from toppling and 
computing all the ways of receiving two grains of sand in a given time step. 

A second source of correction terms is the dissipation rate of sand /. To shed k a , kb grains, 
one must consider all the ways of shedding at least that many, since arbitrarily many grains 
could vanish. These corrections to the q Q d are binomially distributed with parameter /: 



qod(r d a,r d b)= V V q d(k a , k b )(l - f y^ db f k a+ k b -r da -r db f « 

v. - i,- VdaJ 

&a — ' da &b — r db 



h 

Tdb 



since 1 — / is the chance that a grain does not vanish, and r^ a of the grains shed to a and 
Tdb of the grains shed to b must survive. Ideally, the dissipation rate is strong enough that the 
system does not become overloaded with sand, but not too strong that it plays a significant role 
in the avalanche sizes compared to the role that q(0, 0) plays in inhibiting avalanches. Thus, 
as a rule of thumb one should choose / to be one or two orders of magnitude smaller than 
mm{q od (0,0)\o,d G {a, b}}. 

Finally, we define the generating functions {Q d '■ od 6 {aa,ab,ba,bb}} associated to the 
four branching distributions {q a d\od 6 {aa, ab, ba, bb}}: 

oo 

Qod{°da, &db) = ^2 qod{rda,rdb)Vda aa db b for a da, a db € C. 

rda,r db =0 

3.5 Shedding self-consistency equations 

Next we connect the generating functions for the total shedding size V d to the branch distri- 
bution Q d in so-called "self-consistency equations". First, we denote 



a 

V(a) 
Q(a) 
Qod{V{&)) 



[CTaa, Oab, Oba-, &bb), 

{V aa {a),V ab (a),V ba (a),Vbb(v)), 

{Qaa(v), Qab(ff), Qba(v), Qbbtf)), 
QodCPda{a),Vdb^))- 



By the theory of multi-type, multiplicative branching (Galton- Watson) processes (e.g., Eq. 
(10.3) of [25J), we have 

P(a)=a-Q(V(a)), (8) 
where • is the usual dot product. Written more explicitly, 

"Paa — @aaQaa(J~ > aai7~ > ab) i (9) 

Vab = C?abQab(Vba,Vbb), (10) 

V ha = GbaQba{Vaa,V ah ), (11) 

Vbb = VbbQbb{Vba,Vbb), (12) 
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where all the V's are evaluated at a = (a aa , a ab , a ba , Cfefc)- In general, 

'Pod^) = (JodQod(Pda,'Pdb)- 

In words, Eq. ([9]), for example, says that to get an avalanche of sheddings of the four types 
(aa, ab, ba, ba) starting from an aa-shedding, the cascade tree begins with an aa-shedding (that's 
the a aa out front), which causes at the next time step a number of other aa- and a6-sheddings 
distributed according to Q aa , each of which in turn cause numbers of sheddings distributed 
according to V aa and V a b- 



3.6 Two null hypotheses for interacting networks 

To fully understand the effects of interactions between distinct networks (or of modular structure 
within a network) we need to compare the systems described above to an appropriate null- 
model. Below are two different approaches. The first neglects node flavor when calculating 
total avalanche sizes, while the second assumes that the networks are uncoupled. 

Null Hypothesis 1 [Flavorless View] Avalanches in interacting networks should be considered 
in the "flavorless view, " which does not distinguish what portion of the avalanche lies in each 
network. 

The distribution of total avalanche size in the "flavorless view" can be recovered via convolution: 

t 

s* a (t) = J2 s a(i,t-i), (13) 

i=0 

t 

s* b (t) = ^2s b (i,t-i). (14) 



i=0 



Note that here the total avalanche size still depends on in which network the first grain is 
dropped; if it is dropped uniformly at random, then the distribution of avalanche sizes is the 



convex combination of (13), (14) with weights equal to the proportions of a- and 6-nodes. 

The flavorless view of interacting networks is another way of viewing the system; by contrast, 
the "uncoupled null hypothesis" supposes that the networks are not connected at all. 



Null Hypothesis 2 [Uncoupled] The two networks are uncoupled (i.e., not connected by any 
edges). 

Avalanches on uncoupled networks are independent, so the avalanche size distribution on un- 
coupled networks is the product measure of the avalanche size distributions on each network. 

Another null hypothesis for interacting networks, which we do not use here, considers the 
properties of a single network with equivalent size and connectivity (i.e., degree distribution) 
as the system of modular or interdependent networks [32] . 
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3.7 Shedding equations in the flavorless view 



The branching process framework for the flavorless view is straightforward to obtain. Let p a , Pb 
denote the fraction of a- and 6-nodes (i.e., p a = N N ^ Nb , Pb = N -qjjy - ) ■ For o,d G {a, b} let 



Qod(k) :=^2q d(hk - i) 



i=0 

denote the sum of q Q d along the diagonal {(i, k — i) : i E {0, 1, k}}. 

If grains of sand are dropped uniformly at random, then the chance that it lands on an a- 
node is p a (and similarly pb for b- nodes). Thus the avalanche size distribution in the flavorless 
view of the two coupled networks is 

*(*) = PaS* a (t) + p b S* b (t) Vt>0, 

where s a ,sl are the distributions of the total sizes of avalanches begun in a,b, defined in Eqs. 
{13]), ©. 

Consider a grain of sand traveling along an edge; if the graph is fully connected, then the 
flavors of the grain's origin and destination nodes are approximately independent, so we can 
approximate the branching distribution in the flavorless view as 

?(*) = Pll*aa(k) + PaPb{q* ab (k) + + P6<&(*0 V A: > 0. 

These distributions generate 



t=0 k=0 

which are related according to the self-consistency equation 



3.8 Toppling branch distributions u a ,Ub 

Although it is more natural to write down the branch distribution q d of sheddings rather than 
that of topplings, and although the numbers of sheddings tell a different picture about the 
cascade (the total amount of sand shed rather than the numbers of topplings), it is convenient 
to reduce the dimension of avalanche statistics by measuring how many nodes in each network 
toppled in an avalanche rather than how many grains were shed from one network to another. 
As a result, we next derive the corresponding toppling branch distributions u a ,Ub from the 
shedding branch distributions q d, and then we solve for the toppling distributions s a ,Sb- 

The key insight is that a node topples if and only if it sheds at least one grain of sand. Thus 
a grain traveling from a network o to network d topples its destination node with probability 
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1 — g o d(0, 0) . If an a-node topples, what is the chance that it topples t a more a-nodes and t b 
more 6-nodes in the next time step? Denoting this branch distribution by u a (t a ,t b ), we have 

u a (t a ,t b )= Pa(k a ,h)(l- q aa (0,0)) ta q aa (0, 0) k ^ ( °|, 



h a — t a kb — t-b 

k b -t b 



x (l-g a6 (0,0))Sa&(0,0)« 



since the node must have at least t a many a- neighbors and at least t b many b- neighbors, only 
t a ,t b of which topple — which are binomially distributed. Similarly, the branch distribution 
u b (t a ,t b ) is 

u b (t a ,t b ) = Pb(ka,k b )(l-q ba (0,0)) ta q ba (0,0) ka - ta ( t a jx 

x (1-^(0,0)^^(0,0)^^. 

Denoting the associated generating functions by U a {T a , T b ),U b (T a , r b ), we obtain the self-consistency 
relations 

S a = T a U a (S a ,S b ), (15) 
S b = T b U b (S a ,S b ), (16) 

where each S is evaluated at (r a , r b ). 



3.9 Summary of the distributions 

Before solving the self-consistency equations, we pause to summarize the branching process 
framework's distributions and generating functions. (Recall that another common name for 
"branch distribution" is "children distribution" .) 



Table 1: Summary of the distributions and their generating functions 





distribution 


generating function 


degree 


Pa(k a , h),Pb(k a , h) 


G a (u a , U b ) , G b (uj a , LU b ) 


toppling size 


s a {t a , t b ), s b (t a: t b ) 


Sa(T a ,T b ),S b (T a ,T b ) 


toppling branch 


Ua(t a ,t b ),U b (t a ,t b ) 


U a (T a ,T b ),U b {T a ,T b ) 


shedding size 


Pod(r a a,r ab ,r ba ,r bb ) 


'Pod(&aa, Cab, ^bai &bb) 


shedding branch 


qod(rda,rdb) 


Q.od{&da,<?db) 



3.10 Solving the self-consistency equations 



We wish to solve Eqs. (j9|), ((10|), ([llj), (|12j) for V a a,Vab,V ba ,V bb , and Eqs. ([15J), ([16| I for U a , U b . 
Then we obtain their underlying probability distributions by reading the coefficients or by 
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differentiating, since, for example, 

(ta, h) ~ 



1 d ta d tb 
t a W drl a drl b 



=0,7-6=0 



As described in [39 1, since numerical differentiation is prone to machine-precision errors, to 
obtain the best precision it is preferred to use Cauchy's integration formula, 



{ta,h) 



1 



(2mf 



D T a 



dT a dr b 



r. 



integrating over a domain D C C 2 that is the Cartesian product of the largest contours that 
contain the origin and no poles of the generating function. (The generalization of Cauchy's 
integration formula to multiple variables can be found in, for example, Theorem 2.1.1 of |27J.) 

In practice, the self-consistency equations Q-( 12 ), (15 )-( 16 ) are transcendental and difficult 
to invert. However, a generalization of Lagrange's expansion to several variables due to I. J. 
Good [22j provides general tools for inverting self-consistency equations of multi-type branching 
processes. For example, Theorem 9 of [22] allows us to compute coefficients of the distributions, 
one at a time. We state the theorem here for the shedding distributions, but it holds for any 
multi-type branching processes in which each type has a positive chance of giving birth to zero 
children. 



Theorem 1 (Good 1960) Suppose that Q{a) is analytic in a neighborhood of the origin and 
that Qod(0) = Qod(0, 0) 7^ for all o,d £ {a, b}. The probability that the whole avalanche 
consists of exactly m d many od-sheddings, starting from i ol i sheddings, is the coefficient of 



o-. 



m aa -%aa rn ab -i ab m ba -t ba m bb -t bb 



ill! 



'a 



ab 



o-, 



ba 



•a 



bb 



in 



Qm aa Q m ab Q m ba Q m bb 
an. -^nh ~^hn ~~*bb 



-aa -£ a b ~^ba 



o-n dQ, 



where b v u is the Kronecker delta, 



is the determinant, and [i,v run over {aa,ab,ba,bb}. 



(17) 



Even more useful is Theorem 10 of [22J, which explicitly gives the generating function 
solutions rather than just one coefficient at a time. Since we use this theorem to compute 
the toppling size distributions, we state it for the toppling size generating functions. Denote 
t = (t ,t&), S = (S a ,Sb) an d let h(r) = t^t^ 2 for arbitrary nonnegative integers r\,r2- 



Theorem 2 (Good 1960) Under the same conditions as in Theoremyj (namely, U(t) is an- 
alytic in a neighborhood of the origin andU a (0) 7^ 0,^(0) ^ 0), we have 



m a ,m b =0 



m a m b 



m a \m b l 



fim a +m, b 



11 u» 



K=0 



where, as above, fx,u run over the types {a, b}, 5" is the Kronecker delta, and 
determinant. 



(18) 



is the 
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This result holds for any number of types as long as each type has a positive probability 
of being barren; the generalization of (18) to four types (in order to obtain the shedding 
generating functions V d) is straightforward. Taking r\ = l,r2 = in (18) gives <S a , while 
taking r% = 0,r2 = 1 gives S b . In practice, thousands of terms can be obtained by truncating 
the sum in (18) using computer algebra systems. 



4 Examples 

Next we turn to examples, beginning with the easiest (two regular graphs with one-to-one 
coupling) and progressively adding complexity (Bernoulli coupling, other degree distributions) . 
For each example we give the generating functions, and in the next section we compare the 
theoretical predictions of the generating functions obtained from Theorem [2] with numerical 
simulation of sandpiles. 



4.1 Regular (,2 a )-One-to-One- Regular (z b ) 

Arguably the simplest nontrivial intra-degree distribution is the delta function, which yields a 
random regular graph. (Approximating electrical grids by regular graphs is not unreasonable, 
since electrical grids in the United States have been found to have narrow degree distribution, 
namely exponential [H HI [2S])- Similarly, the simplest coupling is "one-to-one": each node 
node in a is connected to exactly one node in 6, chosen uniformly at random. Note that one- 
to-one coupling requires that the number of a-nodes must equal the number of 6-nodes. This 
coupling between networks is hardly realistic for, say, interdependent infrastructure, but the 
degree distributions are so simple — they are products of delta functions — that computing the 
generating functions is easy. Thus one-to-one couplings is a natural first case to solve (e.g., [12J ) , 
but any conclusions for infrastructure or other real systems requires more flexible coupling. 

The degree distributions are 

Pa(ka,h) = S Za (k a )Si(kb), Pb(k a ,h) = 5 X (k a )8 Zb (k b ) , 

which generate 

G a (u] a ,U> b ) = UJa a UJ b , G b (uJ a ,U b ) = aj a U b b - 

The branch distributions simplify to 

qaa(r a a,r ab ) = q b a(r a a,r ab ) = — ^—rS Za (r aa )5i(r ab ), 
q a b(r ba ,r b b) = qbb(r ba ,r b b) = — — r ^(^66)^1^60 ) 

Zb ~r 1 

away from the origin, and 

qaa(0,0) =q ba (0,0) 
q ab (0,0) = q bb (0,0) 



Za + l 

Z b 

z b + l 
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at the origin. The generating functions of the branch distributions are 

z a 1 



Qaa(Vaa,<?ab) = Q ba (cr 

aa ; &ab ) 

Qab(&ba,Pbb) = Qbb{&ba , &bb) 



Zb 1 z b 

Z + 1 2:6+1 



The self-consistency relations Q for the generating functions of the shedding distributions 

Paa, Pab, Pba, Pbb are 



Vaa{3) = (Jaa 

V a b{^) = &ab 

Vba{°) = &ba 

Vbb(<?) = abb 



+ _L_^ (J) , Pat(?) 
_fi_ + _l_^ (?)nt(?) , 



By Theorem [TJ the probability that the whole shedding tree consists of exactly m d many 
od-sheddings, starting from i & many otf-sheddings (for each o, d 6 {a, b}), is the coefficient of 



a aa a ab °ba °^ 



'bb 



111 



Qm aa Qm ab Qm ba Qm bb 
aa ~^ab ~^ba -^-bb 



(19) 



where b v is the Kronecker delta, || • || is the determinant, and [/,, v run over {aa, ab, ba, bb}. For 



two regular graphs with one-to-one coupling, the matrix 5^ — is 
/ 



-l + z c 



^ab^ba 



Za+^aa^ab 



Za+<?laOab 











z b 

_°ab°j£_ 
Zb+&baV b l 
1 

l+z h 
"bb 



-1+X b 



Zb+Cba<7 bb 



z b z b 

Zb+Vbga hb —z b <J ba a bh 



Zb+°baO, 



bb 



Zb+Oba° b 



which has determinant 



-z h a z a a a a ab (-1 + v ba a z b l) + z a (-1 + a z aa a ab ) (-z b + (-1 + z b ) o ha o z h h h ) 



{z a + a'a a a (7 ab ) (z b + O ba O ' 



The product 



m aa0 m abQ m baQ m bb ~ Zb^ab (~1 + Ofeog) + 2 a (-1 + ff^ffqh) (-g, + (-1 + g&) 
^-aa ^ab -^6a -^bb 



2a + CTaaO- afe ) (z& + <7 6a <7^) 
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Thus Expression (19) is 
1 



(1 + Za) (1 + Z b ) 



z a + a aa cr ab \ -l+m aa +n< 



Zb + CTfeaO", 



bb_ \ -l+m ah +m bh 



l + z b 



-ZbO Z a a a (Tab ("I + O-baOfft) + Z a (-1 + 0^0^) (~Z b + (-1 + Z b ) a b a° bb )) ■ 



For arbitrary z a ,z b £ N, we obtain that the toppling branch generating functions are 

{r a + Za) Za {n + z b ) 



U a {T a ,T b ) 



(z a + l)*»(z b + l) ' 
(r a + z a ){T b + z b ) Zb 



(20) 
(21) 



(z a + i)(z b + iy* ' 

Now the toppling size distribution can be read off the generating function obtained from The 
orem[2j which we compute symbolically in Mathematica. 



4. 2 Regular ( z a ) -Bernoulli (p) -Regular (z b ) 



One-to-one coupling is unsatisfactory for describing real world networks, in which nodes in one 
network may connect to arbitrarily many nodes in other networks (even none). A natural next 
choice for inter-network coupling is the Bernoulli distribution: each node has one neighbor in 
the other network (chosen uniformly at random) independently with probability p, no neighbor 
with probability 1 — p. That is, the inter-degrees are Bernoulli distributed with parameter p. 
This allows variable coupling and no longer requires that the number of a-nodes equals the 
number of 6-nodes. The distribution's finite support means that the generating functions have 
finitely many terms. 

For simplicity, we let p be the same for a and b. Then, denoting the probability density 
function of the Bernoulli distribution by B, the degree distributions are 

Pa(k a ,h) = S Za (k a )B(k b ;p), 
Pb{K, h) = B(k a ;p)5 Zb (k b ), 

which generate 

G a (u a ,u b ) = (l-p)u*> +pu* a u) b , 
G b (u a ,uj b ) = (1 -p)u z b h +pu a oj Zb . 

The branch distributions simplify to 
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5 Zb (r bb )5i(r ba ) 



.a(r aa ,r ab ) = —5 Za (r aa )6i(r ab ), 



Za + 1 

q b b(r ba ,r bb ) = 5 Zb (r bb ) 



P 1 Si(r ba ) + - — -5 (r ba ) 
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away from the origin, and 



qaa (0,0) = 1 

q ab {0, 0) 
q ba (0,0) 



P 



1 — p 



Za + l 

Zb 
Z b + V 



Za + l' 



<ftb(0,0) = 1 



■p 



z b + l 



1 — p 

Zb 



at the origin. Note that q ab , q ba are identical to those for Regular-One-to-One-Regular, because 
a grain that is shed between the networks is equally likely to land on any of the nodes that 
have an inter-edge. The generating functions of the branch distributions are 
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The toppling branch generating functions are 

[p(l - T a ) + (z a + \){ Za +T a - 1)] *" (1 + Z b + p{r b - 1)) 



U a {r a ,T b ) 

Ub{T~a,T b ) 



Z Z a a {z a + 1)^(1 + Z b ) 

[ P (i - n ) + (zb + i){z b + n- l)] Zb (i + Z a + p{r a - 1)) 

z z b b (z b + + Z a ) 



When p = 1, these reduce to the generating functions (20), (21) for two regular graphs with 
one-to-one coupling. 



4.3 Current work on degree distributions with infinite support 

Currently we are working on other degree distributions that are not Kronecker delta func- 
tions nor Bernoulli probability density functions. The trouble with Poisson- and power-law- 
distributed degree distributions — and any others that have infinite support — is that the branch 
generating functions contain double-summations that are difficult (perhaps impossible) to sim- 
plify because of the coupling. For example, the shedding branch generating functions for two 
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power-law networks with Poisson-distributed coupling are 



Qaa{kai &&) 



<7aa(0,0) + Yl 



k a +k b >0 




Qab(k a ,h) 



q ab (0,0)+ Yl 



k a +k b >0 



Qba(k a ,h) 



<fta(0,0)+ J] 



fc a +fei,>0 




Qbb(k a ,k b ) 



q b b(0,0)+ Yl 




k a +k b >o 



k a l(k a + k b )((-l + P) 



where £(•) is the Riemann zeta function. Note that the factor l/(k a + kb) couples the double 
summation in a way that prevents us from separating the sums and analytically simplifying 
it. As a result, we cannot simplify the generating function to use a poly logarithm, as you can 
for an isolated power-law network [2T]. Nevertheless, the toppling branch generating functions 
may be analytically tractable, and of course one could always truncate sums before applying 
Theorem [2j 

5 Results of generating function predictions and computer sim- 
ulation 

First we compare the theoretical predictions of the generating function framework to computer 
simulations of the sandpiles. (Another way to match theory and experiment would be to sim- 
ulate percolation with edge traversal probability 1/k, where k is the degree of the destination; 
this type of percolation may resemble processes other than cascading failures in infrastructure.) 
Next we show that the coupling of networks makes them less vulnerable to large avalanches, 
in contrast to the conclusion of the simpler model in [36 1. However, we also show that coupled 
networks suffer avalanches that are large in both networks more frequently than uncoupled 
networks. Finally we report on numerical simulations that are not yet covered by the math- 
ematical framework, in which load is dropped on one network and we measure the size of 
avalanches in the other network, which further illustrates the destabilizing effect of coupling 
between networks. 

5.1 Regular(3)-One-to-One-Regular(10): matching theory and experiment 

Here we compare the theoretical predictions of avalanche size to a simulation with two regular 
graphs a, b with uniform internal degrees z a = 3, Zf, = 10 and one-to-one coupling. The param- 
eters are N a = N b = 10 4 nodes per network, where each node is initialized with a number of 
grains chosen uniformly from to a node's degree minus one (to expedite the simulation) and 
we initially run 10 4 "transient" events (grains for which we do not record statistics). We then 
record statistics for an addition 10 5 grains of sand dropped, using dissipation rate / = 0.001. 
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Figure 1: Comparison of theoretical predictions (red line) with simulation (blue line) of 
avalanche size distributions for Regular(3)-One-to-One-Regular(10). We denote s a (t a ) := 
Ylt^=o s a(t a , tb) (and similarly for s a (tb), Sb(t a ), Sb(tb))- For the theoretical predictions, we com- 
puted 40 terms in Eq. (18), which generate 8000 terms for each generating function. 



The match between simulation and theory is good. As predicted by percolation theory 
(e.g., p2]), the branching process approximation of a random graph of n nodes holds well for 
0(y/n) many nodes. That is, after an avalanche affects approximately 0(y/n) many nodes, the 
"branches" of the "avalanche tree" begin to collide due to the presence of large loops in the 
graph. We see that effect in our simulations: for 10 4 nodes, the avalanche is tree-like up to 
about 10 2 nodes. 



5.2 Coupling between networks can stabilize them individually 

What is interesting about these cascades in interacting networks in Fig. [I] is that, unlike most 
variants of the sandpile model, the avalanche size distributions are not quite power laws. In Fig. 
[2] we plot the generating function predictions with the lines showing the best fits when data 
points 3 through 15 are considered, illustrating that the avalanche size distributions fall short 
of being power laws at large avalanche size. This is not due to finite size effects, because the 
generating function predictions are independent of the number of nodes. Instead the avalanches 
fall short of being a power law for large avalanche size because sand leaks to the other network, 
an effect which is heightened for large avalanches. 
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Figure 2: Avalanche size distributions within the networks are not power laws. Plotted in 
log-log scales in blue and maroon are the generating function predictions of avalanche size 
distributions s a (t a ,tb) (left) and Sb(t a ,tb) (right) for Regular(3)-One-to-One-Regular(10). Plot- 
ted in orange and green are best fit lines to data points 3 through 15, which show that the 
avalanche distributions are sublinear for large avalanches due to sand leaking across to the other 
network (not due to finite size effects). We denote s a (t a ) := ^2^° = q s a (t a ,tb) (and similarly for 
Sa(tb), Sb(t a ), Sb(tb))- Here we computed 40 summands in Eq. (18). 



As illustrated in Fig. [3j when compared to mean-field behavior (avalanche size distribution 
a power law with exponent 3/2), in coupled networks the avalanche size distributions within 
each network have more frequent small avalanches and less frequent large avalanches. Yet, 
the avalanches in the flavorless view (Null Hypothesis [I]) — i.e., the two interacting networks 
viewed as one network without flavors — do show the mean-field behavior, consistent with the 
robustness of the mean-field behavior over different network structures [21J. Nonetheless we 
show that avalanches within a network that is coupled to another are less likely to be large 
than one would predict if the network were isolated. 

This stabilizing effect strengthens with increased coupling between networks. We explore 
this using the generating function predictions for two random regular graphs with Bernoulli- 
distributed coupling (i.e., each node has a neighbor in the other network with probability p). As 
illustrated in Fig. |4j as we strengthen the coupling between the networks (i.e., the parameter 
p of the Bernoulli distribution), large avalanches become progressively less likely while small 
avalanches become more likely. In Fig. [4] networks a and b are random 3-regular graphs 
(internally); we choose the same intra-degree for both networks in order to isolate the effect 
of the coupling, though the same conclusion — that increased coupling between the networks 
mitigates large avalanches and amplifies small ones — holds for random regular graphs with 
different intra-degrees, as well. 

These two effects — that large avalanches are less likely and are further suppressed with 
increasing coupling between networks — suggests that coupling a network to other networks 
stabilizes it, because the other networks act as reservoirs for extra load. This can be understood 
intuitively as follows. Suppose an isolated network is poised to be overwhelmed by an avalanche, 
with many nodes dangerously near their capacities. Coupling the network to another network 
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Figure 3: In networks sparsely connected to others, large avalanches are mitigated and small 
avalanches are amplified compared to the mean-field behavior that the network would exhibit 
if it were isolated (power-law with exponent 3/2), due to the stabilizing effect of connecting 
networks. Here we plot in blue circles and maroon squares the generating function predic- 
tions for the marginalized avalanche distributions s a (t a ), Sb{tb) (defined in Fig. [2js caption) for 
Regular(3)-One-to-One-Regular(10) coupled networks, together with the normalized power-law 
with exponent 3/2 (gold diamonds). 



may mitigate young avalanches since it sheds some of its load during its crucial, early formation, 
before it grows sufficient size and momentum to overwhelm the whole network. Moreover, 
avalanches that leak across to another network are likely to be small there, and only rarely do 
they grow large and hence become likely to amplify the cascade in the original network. This 
stabilizing effect of coupling networks contrasts to the result for a simpler version of this model 
in [36], where they concluded that coupling networks always destabilizes. The assumption in 
|36j that nodes shed sand to every node means that load shed to the other network can only 
exacerbate the cascade; here, load shed to the other network frequently mitigates the cascade, 
and only rarely exacerbates it. 

This result also ameliorates the warnings in [12] about the catastrophic cascades of disrupted 
connectivity in coupled networks. Buldyrev et al. found that coupling two networks exacerbates 
cascades of failing connectivity because coupling the networks can only provide new ways to 
fail and never ways to mitigate them, whereas in this model the coupling can both stabilize and 
destabilize a network. Yet there are regimes where introducing connections between networks 
can destabilize them as shown next. 
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Figure 4: Strengthening the coupling between two networks further mitigates large avalanches. 
Here we show the generating function predictions of the avalanche size distribution s a (t a ) = 
Ylt^=o s a(t a , tb) of avalanches begun in network a for two random regular graphs, each with uni- 
form intra-degree 3 and connected via Bernoulli-distributed coupling with parameter p (shown 
here for p = 0, p = 0.5 and p = 1). As p increases, large avalanches become less likely and 
small avalanches become more likely, due to the stabilizing effect of the coupling to the other 
network. 



5.3 Coupling networks can destabilize them jointly 

Although large avalanches are mitigated in networks that are connected to other networks (com- 
pared to isolated networks), when large avalanches do occur in a network they more frequently 
accompany large avalanches in the other network. The intuition is clear: a large cascade in one 
network likely leaks across the weak coupling to the other network, and the cascades in the two 
networks amplify one another. Said differently, large avalanches in the marginalized avalanche 
size distributions are reduced, but in the joint avalanche size distribution for the two networks, 
avalanches large in both networks become more likely than for isolated networks. 



As in Section 5.2 we compare the joint avalanche size distribution for two coupled networks 
to the null hypothesis of two isolated networks. For two isolated, random regular graphs, 
the avalanche size distributions are mean-field, so their joint avalanche size distribution is the 
product measure of two power-laws with exponent 3/2, 



„ uncoupled 



(t a , tb) 



(tJb) 



-3/2 



(22) 



C(3/2)2 

where £ is the Riemann zeta function. For two coupled networks, the appropriate joint distri- 



27 



bution to compare to Eq. (22) is 



s(t a , h) = g ( s o(*o) *b) + s b (t a , t b )). (23) 



The justification of Eq. (23) is as follows. We drop grains of sand uniformly at random on 
the nodes, and there are as many a-nodes as 6-nodes, so the chance that the sand lands on an 
a-node or 6-node is a fair coin toss (1/2, 1/2). Conditioned on where it lands, the chance that 
the first node topples is roughly the same for both networks, because we find numerically that 
grains are approximately uniformly distributed from to k — 1, where k is the degree of the 
node. Conditioned on these two events, the chance that the ensuing avalanche topples t a many 
a-nodes and % many o-nodes is s a (t a ,tb) or Sb(t a ,tb), respectively. 



Fig. 5 compares the joint distribution (23) of two coupled networks with the joint distri- 



bution (|22j) of two uncoupled networks. The large sea of yellow in the bottom-right corner 
indicates that s(t a ,tb) > s uncoupled (i a , tb) for t a ,tb both large — i.e., that avalanches large in one 
network more frequently accompany avalanches large in the other network in coupled networks 
compared to uncoupled networks. 



5.4 Numerical experiment: dropping grains only in network a 

The multi-type branching approximation above shows that sandpile models on two weakly 
coupled networks rather than on an isolated network reduces the frequency of large avalanches. 
However, dropping sands in only one of the two networks (rather than in both) magnifies the 
distinction between coupled and isolated networks. Here we consider dropping sand only in 
network a, though one could use other rules such as dropping on average two grains in a for 
every one grain dropped in b. 

Dropping sand only in network a pushes a-nodes to their capacities, thereby causing many 
avalanches in a that occasionally leak to b via the sparse coupling. On one hand, network b 
can stabilize network a by serving as a reservoir for dumping extra load. On the other hand, 
avalanches in b can explode in size and thus amplify avalanches in a. Which effect is more 
pronounced asymptotically? 

First we explore how the size and frequency of large avalanches in network b depend on the 
strength of the coupling between the networks. In Fig. [6] we plot the fraction of avalanches 
that topple a certain large fraction of 6-nodes as a function of the coupling strength, which 
here is the mean of the two Poisson inter-degree distributions. There is no phase transition at 
some positive critical coupling strength: once the networks are even slightly coupled (Poisson 
mean 0.01), avalanches in a can overwhelm b. That is, a tiny capacity to leak avalanches across 
the networks suffices to topple nearly all of network b roughly every twentieth avalanche (when 

Z a = 3,Z b = 10). 

Second, we explore how the severity of avalanches in b depend on the relative density of 
intra-edges in a and b. Although the fraction of large avalanches in b saturates once the 
networks are slightly coupled, the size and frequency of large avalanches in b depends on which 
network is more dense (i.e., has more intra-edges). As illustrated in Fig. [6j when a and b 
are random regular graphs with weak Poisson coupling, the large avalanches in b (the network 
not receiving external load) are more frequent and larger when b is less dense (uniform degree 
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Figure 5: Avalanches that are large in both networks (yellow area, bottom-right) are more likely 
in coupled networks than in uncoupled networks. Plotted in color is C[s(t a , h)— s uncoupled (i a , tb)] , 
where C(x) := sgn(x) log(|x| ), and where s(t a ,tb) and s uncoupled (i a , tb) are the joint avalanche 
size distributions for two random 3-regular and 10-regular graphs with one-to-one coupling and 



with no coupling, respectively, defined in Eqs. (23) and (22). Yellow and green correspond 



to s > s uncou P led ; rec j indicates that avalanches of those sizes are equally likely in coupled and 
uncoupled networks; and purple and yellow indicate that s < s uncou P led . Values in the color 
legend are s(t a ,t b ) - s nncou ^ cd (t a ,t b ). 



Zb = 3 < z a = 10, triangles in Fig. [6| compared to when b is more dense (uniform degree 
Zb = 10 > z a = 3, circles in Fig. [6]). 

Intuitively, when a is more dense than b, a has more total capacity — since here capacities 
are degrees — so the large cascades in a consist of so much sand, relative to the total capacity of 
b, that it can easily overwhelm b by jumping across the weak coupling. As a result, to prevent a 
network from large cascades caused by a neighboring network, one should increase the network's 
capacity to be comparable to or larger than the neighboring network's capacity, which in this 
model means add more edges since capacities are the nodes' degrees. 



6 Discussion and Conclusion 

We develop a mathematical framework for interdependent networks and for approximating 
cascades on them using multi-type branching processes. Along the way we elucidate the bias 
induced by requiring matching undirected edge stubs in bipartite and interacting networks. 
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Figure 6: With sand dropped in network a uniformly at random, the avalanches in network 
b are larger when it is less dense than a (triangles) compared to when it is more dense than 
a (circles). We vary the coupling between the networks by tuning the means of the Poisson 
inter-degree distributions from 0.01 to 2.0 (horizontal axis). Plotted vertically is the fraction 
of avalanches that topple > tJVj,, where the threshold t is varied from 0.5 (blue) to 0.95 (gold), 
and Nb is the number of 6-nodes. The networks are regular graphs with N a = Nb = 1000 nodes 
in each network and uniform intra-degree z a = 10, Zb = 3 (triangles, above), z a = 3,Zb = 10 
(circles, below). 



Using a generalization of Lagrange's expansion to several variables [22], we solve the branching 
process equations for the avalanche size distributions for random regular graphs with one-to-one 
or with Bernoulli-distributed coupling, and show the theoretical predictions match well with 
results from simulations. For sandpile models, we find that coupled networks are more stable 
than isolated ones, in that large avalanches occur less frequently due to sand leaking across 
to the other network, in contrast to the conclusion in [36) using a simpler sandpile model and 
to the conclusion in [T2] for a model of cascading failures of connectivity in coupled networks. 
However, we also show that coupling networks can destabilize them, in that large avalanches, 
though more rare, more frequently accompany large avalanches in the neighboring network. 
Furthermore, disparity in capacity and in applied load can enhance large avalanches: when 
load is applied to one network, large avalanches in the second network increase in severity and 
in frequency, an effect that is amplified with increased coupling between the networks and with 
increased disparity in relative capacity. 

These findings suggest economic and game-theoretic implications for infrastructure. On one 
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hand, a greedy owner of an electrical grid (say) prefers to add connections to other networks, 
since this mitigates large cascades in her own network. By contrast, to mitigate avalanches that 
are simultaneously large in the entire collection of connected electrical grids, the grids should be 
less coupled to one another. Thus, like in the game Prisoner's dilemma [43J , the action optimal 
for societal welfare — decrease coupling between networks — conflicts with the action optimal for 
individual networks — increase coupling between networks. More detailed economic models that 
combine results like those here with economic and physical considerations of electrical grids 
may elucidate what networks optimally balance the stability to large avalanches with the cost 
of adding connectivity. 

Here we have focused on mitigating large avalanches, with the example of cascading fail- 
ures in infrastructure in mind. However, cascades in coupled networks apply equally well to 
situations in which we wish to enhance large avalanches. For example, advertisers who design 
word-of-mouth campaigns to spread adoption of products in social networks [8] wish it to spread 
from one sub-population to another across sparse connections. Meanwhile, sociologists who use 
response-driven surveys (RDS) to collect data — in which they pay subjects to recruit friends to 
participate in the survey — wish that their surveys penetrate bottlenecks to spread to different 
populations (say, from drug users to homosexuals) [36]. To understand cascades in interacting 
networks in social contexts such as these would require different branch distributions for the 
branching process (rather than ones inversely proportional to degree, as for the sandpile model 
here) and networks with clustering or with arbitrary subgraphs [29, [3"o] 137] . We expect that 
the tools for solving multiplicative branching processes developed by I. J. Good [22J will find 
other uses for locally tree-like networks, as it is straightforward to implement using computer 
algebra systems; to study cascades on interacting networks that are not tree-like would require 
modification (along the lines of [29 | [35ll37j ). 

Yet there is more work do for the cascades on tree-like graphs considered here, in particular 
to more faithfully capture infrastructure. For instance, different rules for shedding load other 
than "shed one grain of sand to each neighbor" deserve attention (see, e.g., [31]). An example 
shedding rule for interacting networks is that nodes preferentially shed to their network or to 
other networks. Another open problem is to combine models of bearing and shedding load, 
such as sandpiles, with models of cascading failures of connectivity, such as the model in [12] ; 
as we have shown, these two models yield contrasting conclusions regarding the effect that 
coupling between networks has on their stability, so it would be interesting to study which 
effect dominates. 

More broadly, this work suggests that using knowledge about the connection structure 
within and between modules in networks — or alternatively about the structure within and 
between disparate networks — can help to predict processes on thenQ be it cascades, contact 
processes, synchronization, or other dynamics. 
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